library(haven)
library(tidyverse)
library(psych)
library(survey)
library(stargazer)
library(gridExtra)
library(marginaleffects)

schools <- read_dta("DKNS_JBPA_data_final.dta")

## Factor diagnostics

schools <- schools %>% mutate(
  federal_experiment_w2 = -1*federal_experiment_w2 + 5,
  fed_legit_w2 = -1*fed_legit_w2 + 5,
  fed_follow_w2 = -1*fed_follow_w2 + 5,
  fed_vote_w2 = -1*fed_vote_w2 + 5,
  state_experiment_w2 = -1*state_experiment_w2 + 5,
  state_legit_w2 = -1*state_legit_w2 + 5,
  state_punish_w2 = -1*state_punish_w2 + 5,
  state_vote_w2 = -1*state_vote_w2 + 5
)

# Cronbach's alpha

alpha((schools %>% select(federal_experiment_w2, fed_legit_w2, fed_follow_w2, fed_vote_w2)))
alpha((schools %>% select(state_experiment_w2, state_legit_w2, state_vote_w2)), check.keys = T)
alpha((schools %>% select(state_experiment_w2, state_legit_w2, state_punish_w2, state_vote_w2)), check.keys = T)

# Factor scores

fed_unsc <- fa((schools %>% select(federal_experiment_w2, fed_legit_w2, fed_follow_w2, fed_vote_w2)))

state_unsc <- fa((schools %>% select(state_experiment_w2, state_legit_w2, state_punish_w2, state_vote_w2)))

state_nopunish_unsc <- fa((schools %>% select(state_experiment_w2, state_legit_w2, state_vote_w2)))

loadings(fed_unsc)

loadings(state_unsc)

loadings(state_nopunish_unsc)

cor(state_unsc$scores, state_nopunish_unsc$scores, use = "complete.obs", method = "spearman")

schools$fed_unsc <- fed_unsc$scores
schools$state_unsc <- state_unsc$scores
schools$state_nopunish_unsc <- state_nopunish_unsc$scores

# Variable recoding

schools <- schools %>% mutate(
  fed_sc = (fed_unsc - min(schools$fed_unsc, na.rm = T))/(max(schools$fed_unsc, na.rm = T) - min(schools$fed_unsc, na.rm = T)),
  state_sc = (state_unsc - min(schools$state_unsc, na.rm = T))/(max(schools$state_unsc, na.rm = T) - min(schools$state_unsc, na.rm = T)),
  state_nopunish_sc = (state_nopunish_unsc - min(schools$state_nopunish_unsc, na.rm = T))/(max(schools$state_nopunish_unsc, na.rm = T) - min(schools$state_nopunish_unsc, na.rm = T)),
  fed_treatment = ifelse(treat3_new_w2 %in% c(1,3,5,7), "Quality Education", ifelse(treat3_new_w2 %in% c(2,4,6,8), "COVID", NA)),
  prefer_fed = -1*federalism1_w2 + 2, 
  have_child = ifelse(!is.na(housz18), ifelse(housz18 > 0, 1, 0), NA),
  state_treatment = ifelse(!is.na(treat3_new_w2), ceiling(treat3_new_w2/2  - 1), NA),
  state_covid = ifelse(!is.na(state_treatment), ifelse(state_treatment %in% c(1,3), "Covid Claim", "No Covid Claim"), NA),
  state_reason = ifelse(!is.na(state_treatment), ifelse(state_treatment %in% c(2,3), "Knows best", "Lacks resources"), NA),
  ddr = ifelse(ddr_1 < 3, ddr_1*-1 + 3, NA)
)

schools_design <- svydesign(ids = ~1, weights = ~weight_w2, data = schools %>% filter(!is.na(weight_w2)))

# Main results — t-tests, Figure 1

svyttest(state_sc~state_covid, schools_design)

svyby(~state_sc, by = ~state_covid, schools_design, 
      svymean, na.rm = T) %>% ggplot(aes(x = rownames(.), y = state_sc)) + 
  geom_pointrange(aes(ymin = state_sc - 1.96*se, ymax = state_sc + 1.96*se)) +
  ylab("Support for State Gov") + 
  xlab("Treatment") + theme_bw() + ylim(.39, .51)

ggsave("SchoolFigs/StateCOVID.pdf", width = 3.5, height = 5, units = "in")

svyttest(state_sc~state_reason, schools_design)

with(svyby(~state_sc, by = ~state_reason, schools_design, 
           svymean, na.rm = T), 
     state_sc[2] - state_sc[1])/sqrt(svyvar(~state_sc, schools_design, na.rm = T))

svyby(~state_sc, by = ~state_reason, schools_design, 
      svymean, na.rm = T) %>% ggplot(aes(x = rownames(.), y = state_sc)) + 
  geom_pointrange(aes(ymin = state_sc - 1.96*se, ymax = state_sc + 1.96*se)) +
  ylab("") + 
  xlab("Treatment") + theme_bw() + ylim(.39, .51)

ggsave("SchoolFigs/StateReason.pdf", width = 3.5, height = 5, units = "in")

svyby(~state_sc, by = ~state_covid+state_reason, schools_design, 
      svymean, na.rm = T) %>% ggplot(aes(x = rownames(.), y = state_sc)) + 
  geom_pointrange(aes(ymin = state_sc - 1.96*se, ymax = state_sc + 1.96*se)) +
  ylab("") + 
  xlab("Treatment") + 
  scale_x_discrete(labels = c("COVID \nKnows best", "COVID \nLack resources", 
                              "No COVID \nKnows best", "No COVID \nLack resources")) + 
  theme_bw() + ylim(.39, .51)

ggsave("SchoolFigs/State2x2.pdf", width = 7, height = 5, units = "in")

# Appendix B: Table B1 and Figure B1

first_stage <- svyglm(state_sc ~ state_covid*state_reason*fed_treatment, schools_design)

stargazer(first_stage, header = F, covariate.labels = c("State: No COVID Claim", "State: Lack resources", 
                                                        "Federal: Quality education", "No COVID $\\times$ Lack resources",
                                                        "No COVID $\\times$ Quality education", 
                                                        "Lack resources $\\times$ Quality education", 
                                                        "No Covid $\\times$ Lack resources $\\times$ Quality education"), 
          star.cutoffs = c(.05, .01, .001), style = "apsr")

firststage_1 <- predictions(first_stage, 
                            newdata = datagrid(state_reason = c("Lacks resources", "Knows best"), 
                                               fed_treatment = c("Quality Education", "COVID"))) %>% 
  ggplot(aes(x = state_reason, y = estimate, linetype = fed_treatment)) + geom_pointrange(aes(ymin = conf.low, ymax = conf.high), position = position_dodge(.2)) +
  ylab("Support for State Gov") + 
  xlab("Treatment") + 
  scale_linetype_discrete(guide = "none") + theme_bw() + ylim(.38, .52)

ggsave("SchoolFigs/FirstStageReason.pdf", plot = firststage_1, width = 3.5, height = 5, units = "in")

firststage_2 <- predictions(first_stage, 
                            newdata = datagrid(state_covid = c("Covid Claim", "No Covid Claim"),
                                               fed_treatment = c("Quality Education", "COVID"))) %>% 
  ggplot(aes(x = state_covid, y = estimate, linetype = fed_treatment)) + geom_pointrange(aes(ymin = conf.low, ymax = conf.high), position = position_dodge(.2)) +
  ylab("") + 
  xlab("Treatment") + 
  scale_linetype_discrete(guide = "none") + theme_bw() + ylim(.38, .52)

ggsave("SchoolFigs/FirstStageCovid.pdf", plot = firststage_2, width = 3.5, height = 5, units = "in")

firststage_3 <- predictions(first_stage, 
                            newdata = datagrid(state_reason = c("Lacks resources", "Knows best"), 
                                               state_covid = c("Covid Claim", "No Covid Claim"), 
                                               fed_treatment = c("Quality Education", "COVID"))) %>% mutate(label = paste0(state_covid, "\n", state_reason)) %>% ggplot(aes(x = label, y = estimate, linetype = fed_treatment)) + geom_pointrange(aes(ymin = conf.low, ymax = conf.high), position = position_dodge(.2)) +
  ylab("") + 
  xlab("Treatment") + 
  scale_x_discrete(labels = c("COVID \nKnows best", "COVID \nLack resources", 
                              "No COVID \nKnows best", "No COVID \nLack resources")) + 
  scale_linetype_discrete(name = "First Stage Treatment") + theme_bw() +
  theme(legend.position = c(.87,.92)) + ylim(.38, .52)

ggsave("SchoolFigs/FirstStage2x2.pdf", plot = firststage_3, width = 7, height = 5, units = "in")

# Appendix C: Table C1
stargazer(svyglm(state_sc ~ state_covid + state_reason, schools_design), 
          svyglm(state_sc ~ state_covid*state_reason, schools_design),
          svyglm(state_nopunish_sc ~ state_covid + state_reason, schools_design), 
          svyglm(state_nopunish_sc ~ state_covid*state_reason, schools_design), 
          header = F, 
          covariate.labels = c("No COVID Claim", "Lack resources",
                               "No COVID $\\times$ Lack resources"),
          style = "apsr", star.cutoffs = c(.05, .01, .001), type = "html", out = "newc1.html")

# Appendix C: Table C2
stargazer(svyglm(state_experiment_w2 ~ state_covid*state_reason, schools_design), 
          svyglm(state_legit_w2 ~ state_covid*state_reason, schools_design), 
          svyglm(state_punish_w2 ~ state_covid*state_reason, schools_design), 
          svyglm(state_vote_w2 ~ state_covid*state_reason, schools_design),
          header = F, 
          covariate.labels = c("No COVID Claim", "Lack resources",
                               "No COVID $\\times$ Lack resources"),
          style = "apsr", star.cutoffs = c(.05, .01, .001), type = "html", out = "newc2.html")


# Appendix C: Table C3
no_int_state <- svyglm(state_sc ~ state_covid*state_reason + fed_gov_approv_w2 + 
                         state_gov_approv_w2 + cvworry_w2 + 
                         have_child + rol_w2 + l_r_01 + 
                         educ_3_cat + age + female + 
                         ddr + liveRural + 
                         considerAFD, schools_design)

no_int_state_nopunish <- svyglm(state_nopunish_sc ~ state_covid*state_reason + fed_gov_approv_w2 + 
                                  state_gov_approv_w2 + cvworry_w2 + 
                                  have_child + rol_w2 + l_r_01 + 
                                  educ_3_cat + age + female + 
                                  ddr + liveRural + 
                                  considerAFD, schools_design)

stargazer(no_int_state, no_int_state_nopunish, 
          header = F, 
          covariate.labels = c("No COVID Claim", "Lack resources", 
                               "Fed. Gov. Approval", "State Gov. Approval",
                               "COVID Concern", "1+ Children in HH", "Rule of Law",
                               "Conservatism", "Education", "Age", "Woman",
                               "Lived in DDR", "Rurality", "Consider AFD",
                               "No COVID $\\times$ Lack resources"),
          style = "apsr", star.cutoffs = c(.05, .01, .001))

# Appendix C: Table C4
no_int_state <- svyglm(state_sc ~ state_covid*state_reason + fed_gov_approv_w2 + 
                         state_gov_approv_w2 + cvworry_w2 + 
                         have_child + rol_w2 + l_r_01 + 
                         educ_3_cat + age + female + 
                         ddr + liveRural + 
                         considerAFD, schools_design)

no_int_state_supp <- svyglm(state_experiment_w2 ~ state_covid*state_reason + fed_gov_approv_w2 + 
                              state_gov_approv_w2 + cvworry_w2 + 
                              have_child + rol_w2 + l_r_01 + 
                              educ_3_cat + age + female + 
                              ddr + liveRural + 
                              considerAFD, schools_design)

no_int_state_legit <- svyglm(state_legit_w2 ~ state_covid*state_reason + fed_gov_approv_w2 + 
                               state_gov_approv_w2 + cvworry_w2 + 
                               have_child + rol_w2 + l_r_01 + 
                               educ_3_cat + age + female + 
                               ddr + liveRural + 
                               considerAFD, schools_design)

no_int_state_punish <- svyglm(state_punish_w2 ~ state_covid*state_reason + fed_gov_approv_w2 + 
                                state_gov_approv_w2 + cvworry_w2 + 
                                have_child + rol_w2 + l_r_01 + 
                                educ_3_cat + age + female + 
                                ddr + liveRural + 
                                considerAFD, schools_design)

no_int_state_vote <- svyglm(state_vote_w2 ~ state_covid*state_reason + fed_gov_approv_w2 + 
                              state_gov_approv_w2 + cvworry_w2 + 
                              have_child + rol_w2 + l_r_01 + 
                              educ_3_cat + age + female + 
                              ddr + liveRural + 
                              considerAFD, schools_design)

stargazer(no_int_state_supp, 
          no_int_state_legit,
          no_int_state_punish,
          no_int_state_vote, header = F, covariate.labels = c("No COVID Claim", "Lack resources",
                                                                                           "Fed. Gov. Approval", "State Gov. Approval",
                                                                                           "COVID Concern", "1+ Children in HH", "Rule of Law",
                                                                                           "Conservatism", "Education", "Age", "Woman", 
                                                                                           "Lived in DDR", "Rurality", "Consider AFD",
                                                                                           "No COVID $\\times$ Lack resources"), 
          star.cutoffs = c(.05, .01, .001), style = "apsr")


# Appendix D: Balance

with(schools %>% mutate(age_bin = cut(age, breaks=c(10, 20, 30, 40, 50, 60, 70, 80, 90, 100))), 
     chisq.test(state_covid, age_bin))
with(schools, chisq.test(state_covid, ddr))
with(schools, chisq.test(state_covid, l_r_01))
with(schools, chisq.test(state_covid, female))
with(schools, chisq.test(state_covid, liveRural))

with(schools %>% mutate(age_bin = cut(age, breaks=c(10, 20, 30, 40, 50, 60, 70, 80, 90, 100))),
     chisq.test(state_reason, age_bin))
with(schools, chisq.test(state_reason, ddr))
with(schools, chisq.test(state_reason, l_r_01))
with(schools, chisq.test(state_reason, female))
with(schools, chisq.test(state_reason, liveRural))

